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SUMMARY 

Numerical calculations of turbulent reattaching shear layers in a 
divergent channel are presented. The turbulence is described by a 
multiple-time-scale turbulence model. The turbulent flow equations are solved 
by a control-volume based finite difference method. The computational results 
are compared with those obtained using k-e turbulence models and algebraic 
Reynolds stress turbulence models. It is shown that the multiple-time-scale 
turbulence model yields significantly improved computational results than the 
other turbulence models in the region where the turbulence is in a strongly 
inequilibrium state. 


*Work funded under Space Act Agreement C99066G. 
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NOMENCLATURE 

coefficient for tangential velocity correction 
coefficient for transverse velocity correction 
friction coefficient (-r w /( 0. SpU^) ) 
pressure coefficient (*p/(0 . ) 

constant coefficient for eddy viscosity equation (-=0.09) 
turbulence model constants for 6p equation (jG- 1,3) 
turbulence model constants for e t equation (i-1,3) 
constant coefficient (-0.09) 

wall damping function for eddy viscosity equation 
wall damping function for e w equation 
height of backward- facing step 
turbulent kinetic energy (k-kp + k t ) 

turbulent kinetic energy of eddies in production range 
turbulent kinetic energy of eddies in dissipation range 
pressure 

production rate of turbulent kinetic energy 
turbulent Reynolds number (»k^/(i/€^)) 
inlet flow velocity 
time averaged velocity (-{u,v}) 
friction velocity (~%/(r w /p)) 

Reynolds stress (i-1,2,3 and j -1,2,3) 
velocity vector (-{u,v}) 
spatial coordinates (-{x,y,z}) 
reattachment location 
wall coordinate (=u r y/i/) 
deflection angle of the top wall 
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energy transfer rate from production range to 
dissipation range 

dissipation rate of turbulent kinetic energy 

dissipation rate in near-vall equilibrium region 

von Karman constant (-0.41) 

molecular viscosity 

effective viscosity 

turbulent viscosity 

kinematic viscosity of fluid 

turbulent eddy viscosity 

density 

turbulent Prandtl number for kp equation 
turbulent Prandtl number for k t equations 
turbulent Prandtl number for €p equation 
turbulent Prandtl number for e t equation 
wall shearing stress 
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INTRODUCTION 


The experimental study of a reattaching shear layer in a divergent 
channel [1] was designed to test the predictive capability of various 
turbulence models, to identify any deficiency in turbulence closure models, 
and thus to improve predictive capability of turbulence models. The flow 
geometry is shown in Figure 1. The height of the backward- facing step is 
smaller than the boundary layer thickness of the incoming flow. Abrupt 
breakdown of the boundary layer generated a strongly inequilibrium 
turbulent flow. Furthermore, a strong pressure gradient was generated by 
varying the divergence angle of the top wall to study its effect on the 
development of the turbulence field, especially the Reynolds stress, and 
the reattachment process. A number of turbulence models, such as k-e 
turbulence models and algebraic Reynolds stress turbulence models (ARSM) , 
were shown to yield poor computational results for the flow [1,2]. It is 
also shown in References 1 and 2 that a modified ARSM, with modifications 
in the dissipation rate equation, yielded computational results which are 
in good agreement with measured data. However, generality of the improved 
predictive capability for other complex turbulent flows has not been shown 
yet. 

It has been shown previously that the high Reynolds number 
multiple- time -scale turbulence model yields accurate computational results 
for a number of complex turbulent flows such as a wall jet flow, a 
wake-boundary layer interaction flow, a confined coaxial jet without swirl 
and a confined coaxial swirling jet to name a few [3,4]. In the 
single- time-scale turbulence models such as k-e turbulence models, 
algebraic stress turbulence models, and Reynolds stress turbulence models, 
a single time scale is used to express both the turbulent transport and the 
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dissipation of the turbulent kinetic energy. However, this practice is 
inconsistent with physically observed turbulence in the sense that the 
turbulent transport is related to the time scale of energy containing large 
eddies and the dissipation of turbulent kinetic energy is related to the 
time scale of fine scale eddies in the dissipation range. The 
single- time -scale turbulence models yield reasonably accurate computational 
results for simple turbulent flows; however, the predictive capability 
degenerates rapidly as turbulent flows to be solved become more complex. In 
the multiple- time-scale turbulence models [3-7], the turbulent transport of 
mass and momentum is described using the time scale of the large eddies and 
the dissipation rate is described using the time scale of the fine -scale 
eddies. The improved computational results obtained using the 
multiple -time -scale turbulence model for complex turbulent flows can be 
attributed to the physically consistent nature of the turbulence models 
discussed above. 

In numerical calculations of turbulent flows, wall function methods 
are most frequently used to model the near-wall region. These methods have 
been derived from the logarithmic velocity profile based on experimental 
observations that the turbulence in the near-wall region can be described 
in terms of the wall shearing stress. Therefore, these methods are not 
valid if the logarithmic velocity profile no longer prevails in the 
near-wall region. For example, the logarithmic velocity profile no longer 
prevails in the near-wall or in the wake regions of unsteady turbulent 
flows [8], therefore wall function methods can not be applied. Many other 
cases for which the wall function methods are invalid can be found in 
References 9 and 10. Due to this limited applicability of the wall function 
methods, numerous alternative approaches have been proposed. In the 
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alternative approaches, the near-wall low turbulent Reynolds number region 
is included into numerical analyses to overcome the shortcomings of the 
wall functions methods. Various turbulence models which include the 
near-wall low turbulence region can be classified as two-layer (or 
multi-layer) turbulence models [11] and low Reynolds number turbulence 
models [10] based on the way the near-wall region is treated. More detailed 
discussion on the advantages and disadvantages of various near-wall 
turbulence models can be found in References 9 and 10. 

In the present study, the near-wall turbulence is described by a 
"partially low Reynolds number approach." In the model [9], only the 
turbulent kinetic energy equations are extended to include the near-wall 
low turbulence region and the energy transfer rate and the dissipation rate 
inside the near-wall layer are obtained from algebraic equations. The 
algebraic equations were obtained from a k-equation turbulence model [12]. 
It would be appropriate to classify the method as a "partially low Reynolds 
number approach" to distinguish it from other classes of methods. This 
approach was first used in Chen and Patel to solve turbulent flows over 
airfoils [13]. Advantages of the partially low Reynolds number approach 
over the other methods can be summarized as follows. The turbulence length 
scale of the external flows is related to the flow field characteristics 
[14] . On the other hand, the turbulence length scale of boundary layer 
flows is strongly related to the normal distance from the wall. This 
characteristic of the wall bounded turbulent flows can be described quite 
naturally by the present class of turbulence models. The low Reynolds 
number turbulence models can also be used to describe the wall bounded 
turbulent flows; however, more grid points have to be used to resolve the 
steep dissipation rate in the near-wall region. It is also interesting to 
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note that various similar k-equation turbulence models, which form the 
basis of the present near-wall turbulence model, yield accurate 
computational results for a class of simple turbulent boundary layer flows 
[15], turbulent flows with drag reduction [16], and fully developed 
unsteady turbulent pipe flows [8], However, the k-equation turbulence model 
itself is less useful for separated and/or swirling turbulent flows with 
complex geometry due to lack of a systematic method to evaluate the 
turbulence length scale. Development of the near-wall turbulence model and 
its application to fully developed turbulent channel and pipe flows can be 
found in Reference 9. It has been shown in the reference that the present 
near -wall turbulence model can resolve the over -shoot phenomena of the 
turbulent kinetic energy and the dissipation rate in the region very close 
to the wall and that significantly improved computational results for the 
turbulence structure in the near- wall region are obtained. Incorporation of 
the same near-wall turbulence model into a k-c turbulence model and its 
application to complex turbulent flows such as a supersonic turbulent flow 
over a compression ramp and a transonic flow over an axisymmetric curved 
hill can be found in References 17 and 18, respectively. 

The numerical method used herein is based on the pressure correction 
method [19] which has been used most extensively to solve incompressible 
flows the domain of which can be discretized by an orthogonal mesh. 

However, the present numerical method Is applicable for both incompressible 
and compressible flows with arbitrary, complex geometries. The capability 
to solve compressible flows is achieved by including a convective 
incremental pressure term into the pressure correction equation [17,18]. In 
the method, the velocities are located at the same grid points and the 
pressure is located at the centroid of the cell formed by the four adjacent 
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velocity grid points. This grid layout was found to be quite suitable to 
solve flows with complex geometries [17]. The accuracy and the convergence 
nature of the numerical method have been demonstrated by solving a number 
of flow cases. The example problems considered in References 17 and 18 
include: a developing channel flow, a developing pipe flow, a 
two-dimensional laminar flow in a 90 degree bent channel, polar cavity 
flows, a turbulent supersonic flow over a compression ramp, and a shock 
wave - turbulent boundary layer interaction in transonic flow over a curved 
hill. It was found that the numerical method used herein yielded accurate 
computational results even when highly skewed, unequally spaced, curved 
grids were used. 


TURBULENT FLOW EQUATIONS 

The incompressible turbulent flow equations are given as; 
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where eqs. (1-3) follow from the conservation of mass, u-momentum, and 
v-momentum, respectively. In numerical calculation, the conservation of 
mass equation is replaced by a pressure correction equation given as: 
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where the last term represents the mass imbalance, and the first two 
convection terms are unnecessary for incompressible flows. The pressure 
correction equation can be derived following the standard SIMPLE procedure 
[19]. In the present numerical method, all flow variables, except pressure, 
are located at the same grid points and the pressure node has been located 
at the centroid of the cell. The control volume for the pressure correction 
equation is defined as the cell enclosed by the four neighboring grid 
points. The velocity-pressure decoupling is eliminated by treating the 
pressure correction equation as a continuous form partial differential 
equation rather than treating it as a constraint condition. In the former 
case, the discrete pressure correction obtained from eq. (4) becomes a 
five-diagonal system of equations for rectangular grids. On the other hand, 
the discrete pressure correction equation obtained by directly substituting 
the incremental pressure - incremental velocity relations into the 
conservation of mass equation yields a nine-diagonal system of equations. 
The latter discrete pressure correction equation can yield a 
velocity-pressure decoupled solution, whereas the former equation does not 
[17]. 

In control -volume based finite difference methods, the discrete system 
of equations is derived by integrating the governing differential equations 
over the control volume [19]. For curvilinear grids, the number of 
interpolations required to obtain flow variables at the cell boundaries is 
significantly reduced by using the present grid layout. Enhanced 
convergence rate is partly attributed to the grid layout which required 
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fewer interpolations [17]. In solving the discrete system of equations, the 
off-diagonal terms are moved to the load vector term and the resulting 
system of equations can be solved using a tri-diagonal matrix algorithm 
(TDMA) . 


TURBULENCE EQUATIONS 

For clarity, the multiple- time -scale turbulence model supplemented 
with the near -wall turbulence model is summarized below. The turbulent 
kinetic energy and the energy transfer rate equations for the energy 
containing large eddies are given as; 
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where the production rate is given as; 
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The turbulent kinetic energy equation and the dissipation rate equations 
for the fine scale eddies are given as: 
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The turbulent kinetic energy equations, eqs. (5) and (7), are defined for 
the entire flow domain while the energy transfer rate and the dissipation 
rate equations are valid for the flow domain away from the near-wall 
region. The turbulence model constants are given as; a k p-0. 75, a kt -0.75, 
o e p-1.15, a ct -1.15, Cpi“0.21, Cp 2 -1.24, c p 3 * 1.84, c^-0.29, c t2 - 1.28, and 
c t 3*1.66. These turbulence model constants approximately satisfy the 
near-wall equilibrium turbulence condition, the decay rate of the grid 
turbulence [20] , and the turbulence intensity growth rate in a constant 
shear flow [21]. Further discussion on the establishment of these 
turbulence model constants can be found in References 3. 

The energy transfer rate and the dissipation rate inside the near-wall 
layer are given as ; 


where 



C/if V4 k 3/2 
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Note that in eq.(10) represents the standard dissipation rate for 
near-wall turbulent flows in equilibrium state. The dissipation rate given 
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as eq. (9) is formally identical to the one proposed by Wolfshtein [12]. 

For y«0, eq. (9) takes the limit value given as 2i/k/y 2 , which is an 
analytical solution of the turbulent kinetic energy equation for a limiting 
case as y approaches the wall. Slightly away from the wall where the 
turbulence is in the equilibrium state, f c becomes unity. For near-wall 
equilibrium turbulent flows, the production rate (P r ) is approximately 
equal to dissipation rate (c t ) and hence the energy transfer rate (ep) from 
the low wave number production range to the high wave number dissipation 
range has to be approximately equal to the production and dissipation 
rates. Recall that the production rate vanishes on the wall and grows to a 
peak value at y + «15. Hence eq. (9) may not be a good approximation for 
0<y + <15. However, use of the vanishing boundary condition for the 
turbulent kinetic energy on the wall yields a growth rate of turbulent 
kinetic energy and a production rate that are in good agreement with 
experimental data as well as theoretical analysis [9], 

The eddy viscosity away from the near-wall layer is given as; 

k 2 

*t - <Vf- (12) 

€ P 

and that for the near-wall layer is given as; 

k2 

"* C /if f fJL (13) 

n 

where f ^-1 - 1 . /exp ( + A 2 R^ 2 ) is a linear function of the distance from 
the wall in the viscous sublayer and becomes unity in the fully turbulent 
region. A^-0.025 and A 2 -O. 00001 have been used for the near wall layer, see 
Reference 9. The eddy viscosity given as eq. (13) grows in proportion to 
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the cubic power of the distance from the wall. It can be found in Reference 
10 that the near-wall analysis yields the same growth rate of the eddy 
viscosity in the region very close to the wall. However, there also exist a 
few low Reynolds number turbulence models in which the eddy viscosity 
varies in proportion to the fourth power of the distance from the wall, see 
References 9 and 10 for more discussion. 

COMPUTATIONAL RESULTS 

The experimental data for the reattaching shear layer can be found in 
Reference 1. The inlet free stream velocity was 40 m/sec, the boundary 
layer thickness was 0.019 meters, and the height of the backward- facing 
step was 0.0127 meters. The top wall was deflected from -2 degrees to 10 
degrees to generate a strong adverse pressure gradient. 

In numerical calculations, the inlet boundary was located at four 
step-heights upstream of the expansion corner and the exit boundary was 
located at approximately 35 step -heights downstream of the expansion 
corner. The flow domain was discretized by a 105 by 85 mesh with 
concentration of grid points near the expansion corner and in the bottom 
wall region, see Figure 2- (a). The grid in the vicinity of the expansion 
corner is fine enough to resolve details of the large eddies subjected to 
strong shear and sudden expansion, see Figure 2-(b). The inlet boundary 
conditions for the tangential velocity, the turbulent kinetic energies, and 
the dissipation rates (Cp and € t ) were obtained from experimental data for 
a fully developed boundary layer flow over a flat plate [3,22], The 
non-dimensional velocity and the turbulent kinetic energy profiles were 
scaled to yield a boundary layer thickness of 0.019 meters at the inlet 
boundary. The no-slip boundary condition for velocities and vanishing 
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turbulent kinetic energies were prescribed at the solid wall boundary. At 
the exit boundary, a vanishing gradient boundary condition was used for all 
flow variables except the pressure. A uniform pressure was prescribed at 
the exit boundary. The partition between the near-wall layer and the 
external region was located at approximately y + -100, 12 grid points were 
allocated inside the near-wall layer. The mesh size of the first grid point 
on the bottom wall was Ay + ^2 and the grid size in the normal to the wall 
direction was increased by a factor of approximately 1.15. Further details 
on the computational procedure can be found in Reference 17. 

The calculated streamline contours are shown in Figure 3. The flow 
field consisted of two recirculation zones. The primary recirculation zone 
extended from the separation corner toward the downstream direction; and 
the secondary recirculation zone was very small and confined in the corner 
region. The reattachment location versus the top wall deflection angle is 
shown in Figure 4. It can be seen in the figure that the k-c and ARSM 
turbulence models largely under-predict the reattachment location. The 
modified ARSM yielded a significantly improved computational result, 
however, the present computational result compared more favorably with the 
measured data than did the modified ARSM. 

The static pressure contour lines are shown in Figure 5, where the 
pressure has been normalized by the inlet total pressure and the 
incremental pressure between the contour lines is 0.005. It can be seen in 
the figure that a few contour lines pass through the expansion corner, and 
thus there exists a mild base pressure in the backward- facing step region. 
The calculated static pressure on the wall is compared with experimental 
data as well as the numerical results of Reference 1 in Figure 6. The mild 
pressure drop at x/H-*0 represent the base pressure. For a-0®, the present 
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computational result compared slightly more favorably with the measured 
wall pressure than the other computational results. For a-6°, all the 
computational results compared decently with the measured data; however, 
the slope of the wall pressure in the continuously diverging downstream 
region obtained in the present study compared more favorably with the 
measured data than the other results. This difference may due to the 
different numerical methods used. 

The calculated wall shearing stresses are shown in Figure 7. It can be 
seen in the figure that the location of the peak wall shearing stress 
obtained using the k-e turbulence model is grossly in error. It is 
interesting to note that the modified ARSM under-predicts the peak value 
and the present turbulence model over-predicts the peak value even though 
the relative differences are almost the same for both deflection angles. 

The mean velocity, the turbulent kinetic energy, and the Reynolds 
stress profiles at four downstream locations are compared with experimental 
data and with the calculated results using the modified ARSM [1,2] in 
Figures 8-10, respectively. The experimental turbulent kinetic energy shown 
in Figure 9 was estimated using the measured value of u # ^+v'^ an d an 
assumption that w' ^-(u' 2+ V ' 2 )/2 . As shown in Figures 8-10, both 
computational results exhibit fair comparison with the experimental data. 

It can be seen in Figure 9 that the peak value of the turbulent kinetic 
energy and the shape of the turbulent kinetic energy profile obtained using 
the multiple- time -scale turbulence model compare slightly better with the 
measured data than those obtained using the modified ARSM at x/H-1.0 where 
the turbulence is in a strongly inequi librium state. It has been shown 
previously that the improved computational results for complex turbulent 
flows are attributed to the capability of the multiple- time-scale 
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turbulence model to resolve the inequilibrium turbulence [3], The same 
argument can be applied for the present flow case. At further downstream 
locations, the present computational results compared slightly less 
favorably with the experimental data, a result perhaps due to the near-wall 
turbulence model which can not take into account of the inequi librium 
turbulence . 

The ratio of turbulent viscosity to molecular viscosity at three 
downstream locations are shown in Figure 11. It can be seen in the figure 
that the Jones-Launder k-« turbulence model over-estimates the ratio so 
that the reattachment location is largely under-predicted. On the other 
hand, the present computational results compare quite favorably with the 
measured data so that the reattachment location is correctly predicted. The 
calculated production and dissipation rates of the turbulent kinetic energy 
at the same downstream locations were qualitatively and quantitatively 
almost the same as those of reference 1. 

CONCLUSIONS 

Numerical calculations of reattaching shear layers in a diverging 
channel using a multiple- time -scale turbulence model supplemented with a 
near-wall turbulence model have been presented. The calculated reattachment 
location versus the top wall deflection angle obtained using the present 
turbulence model was in excellent agreement with measured data. The 
calculated wall pressure and the wall shearing stress were also in good 
agreement with the experimental data. The rest of the present computational 
results such as the normalized velocity profiles and the Reynolds stress 
profiles compared favorably with experimental data. The computational 
results obtained using the multiple- time -scale turbulence model compared 
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slightly more favorably with the experimental data than those obtained 
using the modified algebraic Reynolds stress turbulence model. It has been 
shown that prediction of the correct reattachment location depends on the 
prediction of the correct level of the turbulent viscosity, which depends, 
in turn, on the capability of a turbulence model to resolve the entire 
turbulence structure of the flow field correctly. Thus the improved 
computational results are attributed to the capability of the 
multiple- time-scale turbulence model to resolve the strong inequilibrium 
turbulence in the vicinity of the expansion corner and in the following 
shear- layer region. 
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